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A technique to study collisionless dynamics of a homogeneous superconducting 
system is developed, which is based on Riccati parametrization of Wigner distri- 
bution function. The quantum evolution of the superconductiung order parameter, 
initially deviated from the equilibrium value, is calculated using this technique. The 
effect of a time-dependent BCS paring interaction on the dynamics of the order 
parameter is also studied. 

I. INTRODUCTION 

In this work we study the dynamics of the superconducting order parameter within the 
Wigner distribution function approach. The problem of nonstationary phenomena in super- 
conductors has been attracting attention for a long time^*^. The general method for descrip- 
tion of nonstationary and nonequilibrium processes is the Keldysh technique for nonequi- 
librium real time Green's functions^. The equations for superconducting Keldysh Green's 
functions^ are a set of quite complicated nonlinear integro-differential equations, which are 
nonlocal in time and space domains. These equations are considerably simplified in the 
quasiclassical approximation by integrating Green's functions over = /2m — /x,® The 
quasiclassical Larkin-Ovchinnikov equations are still nonlocal in time, but are local in space. 
In the stationary case, these equations transform into Eilenberger's equations^, which are 
effective tools for solving stationary inhomogeneous problems. 

When the time-dependent processes in superconductors are considered, three time scales 
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are the most essential. The time tp ~ u!~^ {up is a plasma frequency) characterizes the 
scale at which the self-consistent scheme for the electromagnetic fields A(r, t), 0(r, t), and 
for the BCS pairing field A(r, t) is established. The time to ~ is the superconduct- 

ing gap) is an intrinsic time for superconductors, during which quasiparticles with energy 
spectrum y^AM-^ are formed in the superconductor. The stage of the relaxation of a 
nonequilibrium disturbance in the quasiparticle distribution is determined by the energy 
relaxation time r^, which is caused by electron-phonon inelastic processes. For conventional 
superconductors, at temperature T, not too close to the critical temperature Tc, we have a 
hierarchy of the characteristic times: tp <^ to <^ t^. At the time interval t ~ ^ to the 
superconductor's dynamics is described by the quasiclassical Boltzman kinetic equation for 
the quasiparticle distribution function together with a self-consistent equation for A{r,t) 
(Aronov-Gurevich equations^). In the opposite case t <^ r^, the dynamics of the supercon- 
ducting order parameter should be described by the quantum kinetic equation. Considering 
the coUisionless evolution of the superconducting order parameter (t <ti r^), the equations 
for the Keldysh Green's functions are reduced to simpler equations for the Green's functions 
at coinciding times. The latter can be transformed to the quantum kinetic equation for the 
Wigner distribution function (WDF). The coUisionless kinetic equation for superconducting 
WDF can also be obtained starting directly from the nonstationary generalized Hartree-Fock 
approximation for the BCS pairing modelS^ (see alsoiSiii). 

Wigner introduced a distribution function in the phase spaced as a quantum analog of 
the classical Boltzman distributions. To study quantum transport, the Wigner- function for- 
malism possesses many advantages. It is extensively used for the description of normal metal 
and semiconducting electron devices whose behaviour is dominated by quantum interference 
effects, e.g. for self-consistent treatment of transient response to a change in the applied 
voltag In recent years, Wigner functions are widely used in the field of quantum optics 
to describe the effects of quantum coherence and superposition in optical systemsti^. Such 
effects are of great interest in qubit (quantum bit for quantum computation) investigations^^. 

CoUisionless dynamics of the superconducting order parameter has gained renewed at- 
tention after the discovery of the BCS-like paired state in dilute fermionic gases^^. The 
ability to control and change the strength of the pairing interaction in these systems opens 
possibilities for new experimental investigations of the dynamics of the order parameter. Re- 
cently, time-dependent BCS pairing was studied theoretically in Ref.-^. The WDF technique 
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developed in this paper provides a useful tool for studying such problems. 

In Sec. II, following to Kulik's approach^, we derive a quantum kinetic equation for 
superconducting WDF in (r, t)-space. This equation is simplified for the case of homogeneous 
state (Sec. Ill) and then used to study the collisionless dynamic of the order parameter in 
small superconducting systems (Sec. IV). The problem of the time evolution of the order 
parameter, initially deviated from the equilibrium value, is considered. It appears that on 
times much smaller than r^, the time dependence of A has an oscillatory nature. Earlier, 
such problem was studied by other authors using a linear response approachi^, assuming 
small deviation from equilibrium. In the present work, these dependencies were obtained 
under arbitrary initial conditions (not only small). The time dependent response of the 
order parameter to a time- varying pairing potential is also studied. A numerical method for 
solving equation for WDF, which is based on Maki-Schopohl transformation^^, is developed. 



II. WIGNER DISTRIBUTION FUNCTION FORMALISM FOR THE 

SUPERCONDUCTING STATE 

We write the Hamiltonian of the superconductor as 

H = Hq + Hi, 

where Hq includes electron interactions with external fields, the vector potential A(r) and 
the scalar potential f{r), as well as with the pairing field A(r), 



Ho=Y. ^r^iW [e-/x + e^(r)]^.(r)- /dr [A(r)^j(r)^j(r) + A*(r)^i(r)^t(r) 

2 



1 

2m 



:A(r) 



(2) 



(we use the system of units where h = ks = 1). Here V'o-li") = S ^pcr(^)6*^'" is the 

p 

annihilation operator of an electron with the spin cr, and fi is the chemical potential. The 
Hamiltonian Hi describes impurity, electron-phonons, electron-electron etc. scattering pro- 
viding processes for relaxation. 

The pairing field A(r) is to be determined from the self-consistency equation 



A*(r) = yo(^|(r)^|(r 



(3) 
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where Vq is the pairing potentiaL The electromagnetic potentials obey Maxwell equations, 

V X A(r) = ^j(r), (4) 
1 d 

VV + -^V- A = -47rp(r), (5) 
where p(r) and j(r) are the charge and current densities, respectively: 

pir) = eJ2{i^tir)M^)), (6) 



.2 



= — E i^imMr) - {WAr))Mr)) - — A(r) ^ {^Ur)Mr)) ■ (7) 
m ^ — ^ mc ^ — ^ 

By introducing the "particle-hole" (Gor'kov-Nambu) representation of the electron cre- 
ation and annihilation operators in terms of 2- vectors 

^p=f V ^l={al,a^,A. (8) 



v&(r)= ( ), ^\r) = Ui^r) ^,(r) ) , (9) 



we define the matrix /pq(t) in the "particle- hole" space. 



where angle brackets denote statistical averaging and = 1,2 are the indices of the 
vectors Ap. The function f^^ is the Fourier transform of the Wigner distribution function 
(WDF) /q,^(p, r, t) generalized to the superconducting case. 



Mp, r, t) = E e^^"" {Al_^Jt)Ap_,^Jt)) . (10) 
q 

Correspondingly, the components of the matrix /(p, r, t) are expressed through the Nambu 
operators \1/Q(r, t) in the Heisenberg representation as 

U = I rfr'e-P^'(v[/t(r + r72,t)vI/^(r-rV2,t)>. (11) 

It follows from Eq. (fTTj) that /u and are real functions, and /12 = /2i- The self-consis- 
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tency equations, Eqs. (jSl), (©, and ((Zj) can be written in terms of / as 

dp 



A = Vo 
p = e 
j = 



(27r)3 
dp 



Trr_/(p) 



(27r)3 
dp 



Trr3/(p), 



i2ny 



Trp/(p), 



(12) 
(13) 
(14) 



where p = p — er^A/c, r_ = (1/2) (ri — 2x2), and Tj are the Pauh matrices. 

The evolution equation for the WDF can be derived from the equation of motion for the 
electron field operators ip = ilJair^t): 



(15) 



Restricting our consideration by the coUisionless stage of the evolution, we neglect the in- 
teraction part Hi of the Hamiltonian and obtain from Eq. (fT3j) the equations of motion for 
the Nambu operators ^E'(r, t) 



^ = 0, A 




(16) 



where ^ = — (V+ier3A/c)^/2m — /i. By making use of the definition of the WDF in Eq. (Illj). 
we arrive, after some algebra, at the following dynamic equation for /(p,r,t) 

(p-^V/2)2 



df ^. 



2m 



+ 1 



eipT3 - A, / 



0, 



(17) 



where [. . .] denotes usual commutator, in which we consider V as an integral operator with 
the kernel Vr5(r — r'), thus (V/) = —(/V) = V/. The quantity [. . .j® is defined as 
[y4, 5]^ = B — B ^ A, where {A ® B) (p, r, t) is the Fourier transformation of the spatial 
convolution {AB){ri,r2) = J dr A{ri,r)B{r,r2): 



(A® 



B){p, v) = j dr'e~'P'''{AB){r + r' /2, r - r'/2) = eH^^^'^-^P '^r^l^Kp, r)5(p, r). (18) 



By making use of the transformation / exp(ir3x/2)/ exp(— 'ir3x/2), we can exclude the 
phase X of the superconducting order parameter and proceed to gauge invariant quantities, 
i.e. the momentum of the superfluid condensate p^ and the potential $ defined by 



^ 2e ^ 
mVs = - Vx A 



(19) 
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The electromagnetic fields are related to p<j and $ through 



eE 



dp, 



's 



- V$ 



eH 



-V X p^. 



(20) 



dt 



This results in the substitutions p p + r-^Ps and ev? — ^ $ in the dynamical equation (|T7jl . 
as well as in the definition of the current in Eq. (fT^ . Note that the anisotropic term p ■ 
arising from p in Eq. p7|) commutes with / and thus drops out from this equation. 

While the physical sense of Ps is obvious, the interpretation of the gauge invariant po- 
tential $ is less evident. Within the framework of the two-fiuid model, it can be interpreted 
as the difference $ = /i^ — /^n between the electrochemical potentials of the condensate of 
Cooper pairs, fig = fi + {l/2)dx/dt and quasiparticles fin = fJ^ — ^ip; thus, a nonzero value of $ 
means the nonequilibrium of the electrons in the superconductor. In bulk superconductors, 
$ and Ps decay within their corresponding lengths: London (skin) depth 6 for p^ and the 
electric field penetration depth A^; for $. 

III. WIGNER DISTRIBUTION FUNCTION FOR HOMOGENEOUS 
SUPERCONDUCTING SYSTEMS 

In what follows, we focus on homogeneous superconducting systems in pure limit, as- 
suming the scattering rate is much smaller than A. To be more specific, we assume the 
magnitude of the order parameter A and the gradient of its phase, Vx, to be uniform inside 
the superconductor. Using an appropriate gauge transformation, we include the spatially 
varying part of the phase of A into the homogeneous p^. A "residual" spatially uniform 
phase is kept to describe the dynamics of the phase of the order parameter. It can be related 
to, e.g., possible (time-dependent) phase on either sides of a Josephson junction. In this 
case, the equation for the WDF takes the form. 



where C,p = + ^ + mVg/2 and C,p = /2m — fi. The phenomeno logical collision term 
z/(/ — /o) qualitatively describes slow relaxation of the WDF to its equilibrium value /o 
which is associated with the interaction Hamiltonian Hi. In the coUisionless limit considered 
below, we will assume u — +0, in order to provide correct analytical behavior of the WDF 
at t ^ +00. 



dl 
dt 



+ z[epr3-A,/] + z/(/-/o) = 0. 



(21) 
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Equation ()2ip has several important properties which can be derived from the equations 
for the matrix elements of /, 

.dfn 



' dt 

. Of 12 

' dt 

. Of 21 

' dt 



2^/12 + A(/n-/22), 
24/21 + A*(/n-/22). 



(22) 
(23) 
(24) 



First, we note that only the difference /n — /22 of the diagonal components of the matrix 
/, enters the equations for the off-diagonal components fu and /21. Furthermore, from 
Eq. (j221), one finds that the sum of the diagonal components fu + f22 = const. This allows 
us to present the function / in the following form: 



/ 



2 L 



l(l-^_)-/X 



/ 




(25) 



where / and g are isotropic functions, and the time-independent quantities J-'± have the 
meaning of quasiparticle distribution functions which are conserved during the stage of the 
collisionless evolution. Assuming the system to be initially in equilibrium and comparing 
Eq. (j25j) with the equilibrium form of the WDF, which can be directly obtained from the 
definition in Eq. (jlOj) 



7o = ^ { 1 ( 1 - ) - ^ ( e^ra - A ) 



(26) 



we find the distribution functions 
1 



± 



tanh^^^±P^±tanh^ 



P-v,(0) 



2T 



2T 



ll + |AP, 



(27) 



and the equilibrium values of the functions / and g 



fo 



A 



9o 



(28) 



In this representation, the dynamic equation ()2H) for the WDF reduces to the following 
system of scalar equations for the functions g and /: 



dt 



2i{Ag-i,f), 



(29) 
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which, together with Eq. ()28p . lead to the normahzation condition 

+ //* = 1. (30) 
The self-consistency equation has the form 

A(t)=i y 177^ rfep/(ep,t)n, (31) 

where ud is the Debye frequency, A = A^(0)Vo is the dimensionless pairing constant, A^(0) 
is the electron density of states per spin at the Fermi level, and Qp denotes angle variables 
associated with the momentum vector. The charge and current densities are given by 

Pit) = -eiV(O) d^pgi^p,t)J^+, (32) 

-^pjd^pj^_, (33) 

where n is the net electron density. Equation ()33|) shows that the electric current is governed 
directly by the superfluid velocity and has nothing to do with the evolution of the WDF, 

j(t) = env,(t) + e(n, - n)v,(0) = j(0) + en[v,(t) - v,(0)], (34) 

where Ug is the condensate density calculated for the initial superfluid velocity Vs(0). This 
property reflects the specifics of the collisionless regime, in which the normal component 
of the current flow is not affected by scattering, and therefore the velocities of both the 
superfluid and normal components of the electron fluid undergo equal changes Vs(t) — Vs(0): 
Vs(0) Vs(t), v„(0) = ^ Vs(t) — Vs(0). From this we conclude that at nonzero tem- 
perature, when the density of the normal component, n„ = n — n<j, is nonzero, the current 
reverses its direction with respect to Vs(t) if the latter becomes smaller than Vs(0)n„/n. 

IV. COLLISIONLESS EVOLUTION OF THE ORDER PARAMETER IN 

SUPERCONDUCTORS. 

In the paper by Volkov and Kogan^^, the problem of evolution of the order parameter 
A(t) at given initial value of the WDF (and corresponding initial self-consistent value of 
A = A(0)) was analyzed within a linear approximation, assuming small deviations of A(t) 
and /(^,t) from their equilibrium values . It was shown that the time variations of A have 
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the form of harmonic oscillations with the period of the order of and the amplitude 
slowly decreasing as t~^/^. At large t ^ = A^^(O), the order parameter approaches a 
constant value Aqo = A{t —>■ oo), which is determined by the initial conditions and coincides 
neither with A(0), nor with the equilibrium value Aq. 

In this paper, we address a more general nonlinear problem, with arbitrary initial con- 
ditions, which may essentially differ from the equilibrium state. In particular, this allows 
us to consider formation of the superconducting state from the initial normal state at low 
enough temperatures, or destruction of the initial superconducting state at high tempera- 
tures. To this end, we apply a numerical procedure, by making use of the so-called Riccati 
parametrization of the functions / and g. Due to the normalization condition (jHOI), these 
functions can be expressed through a single function a{^p,t), 

9 = / = (35) 

1 + aa* 1 + aa* 

which satisfies a nonlinear Riccati-type equation, 

^ = z (-24a - A*a' + a) . (36) 
In the stationary limit (A = const), the solution of Eq. (j36p is 

«o = J——- (37) 

Sp + 



In a general non- stationary case, one needs to integrate Eq. fl36|) . together with the 
self-consistency equation ()31|). Thus, proceeding to the discrete time variable, t = nSt, 
n = 0, 1, . . ., one has to calculate the new value of A from Eq. (jHT|) after each time step 
6t, and then use it for the next step. For sufficiently small 6t, A can be approximately 
considered as constant between t and t + 6t, which allows us to apply an analytical solution 
of Eq. (|36|) within this time interval, 

A*(t)a(t) + — iep cot{ep6t) 
and thus to calculate a(t + 6t) explicitly. As the result, the numerical procedure reduces to 
the numerical solution of the self-consistency equation at each step of calculations. 

In our calculations, we use time steps 6t = 0.02to- After each step, the values of the 
modulus and the phase of A(t) were recalculated by means of the self-consistency equation 
(31). In Figs. 1 and 2, we present time variations of the order parameter modulus with 
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FIG. 1: Collisionless time evolution of the order parameter with initial value A(0) larger than 
the equilibrium value Aq at T=0. In all figures the time is normalized on to = I/Aq- 




5 10 15 20 25 30 



FIG. 2: Collisionless time evolution of the order parameter with initial value A(0) < Aq at T = 0. 

the initial values A(0) essentially different from the equilibrium value Aq at T = 0. It is 
obvious that equal values of A(0) may be obtained for different forms of the initial Wigner 
distribution function /(O). In our evaluation, we use the equilibrium form of /(O) given by 
Eq. fl2(ij) at T = 0, with a formal parameter Ajn, which, however, appears to be slightly 
different from the initial self-consistent value A(0). This difference weakly depends on the 
value of the pairing constant A, which we in the following put A = 0.5. The initial value 
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of Ain = 1.5Ao leads to A(0) ^ I.SAq (Fig. 1), whereas Ain = O.SAq yields self-consistent 
A(0) ^ 0.67Ao (Fig. 2). 

Another type of perturbation in the system is the switching of the A from one value to 

another. Or, more generally, the case of time dependent BCS pairing. We have used the 
equations (35), (36), (38) and (31) (with A = A(t)) to study this problem numerically. In 
Fig. 3 the coUisionless evolution of the order parameter under the changing of A in time is 
shown. 

3.6 
3 
2.5 
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5 10 15 20 25 30 35 40 45 50 
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FIG. 3: CoUisionless time evolution of the order parameter under the changing of coupling constant 
A from the value 0.5 to the 1. 

It is interesting to note that the initial BCS form of the WDF automatically leads to 
conservation of arbitrary initial value of the order parameter phase Actually, such prop- 
erty is associated with the definite symmetry of the initial WDF with respect to ^p, which 
holds during the time evolution, f{^p,t) = f{—^p,t), g{^p,t) = — (?(— ■^pjt), and manifests 
equality of the populations of the electron- and hole-like excitations with equal energies €p. 
The introduction of an imbalance between the electron and hole branches of the excitation 
spectrum (i.e., violation of the above-mentioned symmetry) produces an excess charge in 
the quasiparticle subsystem which, due to electroneutrality of the metal, should be com- 
pensated by the opposite charge of the superfluid condensate. This means the appearance 
of the difference 6fi between the electrochemical potentials fin and fis of excitations and 
the condensate, respectively, which produces time variations of the order parameter phase 
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according to the relationship dx/dt = 2Sfi. For a given constant Sfi, we find continuous 
variation of the phase with a constant rate. 

The processes of formation and destruction of the superconducting state can be also 
analyzed within the nonlinear collisionless approach. By starting evaluations from a very 
small value of Aj„, (~ 10"'^ Aq) in Eq. (j^H|l at T = 0, which approximately represents initial 
normal state (see Fig. 4, we observe a rapid increase in A{t) at the time t ~ to up to 
A ~ Aq, followed by an oscillatory approach to a stable superconducting state. We note 
that the asymptotic value A^o appears to be noticeably lower than Aq, which means that the 
real equilibrium value of A at the superconducting transition is formed via the relaxation 
processes. 

1.5 1 , , , , , 1 
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FIG. 4: Instability of the equilibrium normal state at T = 0. We start from A(0) = O.OOIAq. 

Strictly speaking, at any temperature, including the region T < Tc, the self-consistency 
equation Q always has a trivial solution A = 0, which corresponds to the normal state. 
However, at T < T^, the normal state is associated with a maximum of the free energy, and 
therefore Fig. 4 actually illustrates the thermodynamic instability of the normal state with 
respect to an infinitesimal A, which develops through the quantum kinetic process described 
by Eqs. (j^^ and 0. It is interesting to note that, despite the strong nonlinearity of the 
process, the oscillations of A(t) have almost purely harmonic shape. 

The instability of the superconducting state at the temperature T > Tc is illustrated by 
Fig. 5, which was obtained by starting evaluations from the initial superconducting state in 
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Eq. ()28|) at high enough temperature T = 2.5Ao. The order parameter decreases approxi- 
mately exponentially with the characteristic decay time 0.42to without any oscillations. At 
the final stage of the evolution, the order parameter enters the fiuctuation regime which is 
out of the framework of our self-consistent approach. 

0.4 
0.3 

A/Ag 0.2 
0.1 





0.2 0.4 , 0.6 0.8 1 

FIG. 5: Instability of the superconducting state at T = 2.5Ao > Tc. Initial condition A(0) = 
O.SlAo. 

In conclusion, the authors are grateful to B.Z. Spivak and A.M. Zagoskin for discussing 
problems encountered in this work. 
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